These demonstrations go beyond the content of the course. I don't expect you to understand what is going on in these snippets of code. But I show them to you so you see what sorts of things can be implemented by manipulating the Fourier coefficients.
import numpy as np
from numpy.fft import fft, ifft, fft2, ifft2, fftshift, ifftshift
import matplotlib.pyplot as plt
%matplotlib inline
plt.ion()
from scipy import ndimage, misc
from scipy.signal import gaussian
from IPython.display import display, HTML
display(HTML("<style>.container { width:98% !important; }</style>"))
np.set_printoptions(edgeitems=30, linewidth=250,
formatter=dict(float=lambda x: "% .4f" % x))
def Gaussian2D(dims, std):
'''
Creates an array containing a 2-D Gaussian function, with mean
at mu, and standard deviation std.
Inputs:
dims is a 2-vector containing the number of rows and columns
std is the standard deviation
'''
nr, nc = dims
rr, cc = np.mgrid[0:nr,0:nc]
ctr = GetCentre(rr)
temp1 = rr - ctr[0]
temp2 = cc - ctr[1]
temp1 = np.exp(-temp1**2/std)
temp2 = np.exp(-temp2**2/std)
blah = temp1*temp2
return blah/sum(blah.flatten())
def GetCentre(f):
return np.array(np.floor(np.array(np.shape(f))/2.), dtype=int)
# If f is real-valued...
f = np.random.random(8)
print(f)
[ 0.6558 0.3726 0.1559 0.8886 0.5927 0.6509 0.9460 0.7658]
# Fourier coefs will be cojugate symmetric
F = np.fft.fft(f)
print(F)
print(F[2], F[-2])
[ 5.02823474+0.j -0.22050572+0.90007346j 0.1466144 +0.63084762j 0.34663392-0.68017999j -0.32749292+0.j 0.34663392+0.68017999j 0.1466144 -0.63084762j -0.22050572-0.90007346j] (0.146614397571887+0.6308476194085022j) (0.146614397571887-0.6308476194085022j)
rr, cc = np.mgrid[-128:128, -128:128]
%matplotlib inline
d = np.sqrt(rr**2 + cc**2)
circle = np.zeros([256,256])
circle_mask = d<15
circle[circle_mask] = 1.
circle = ifftshift(circle)
circle[0:22,0:5] = 1.
plt.figure(figsize=(8,4))
plt.subplot(1,2,1); plt.imshow(fftshift(circle), cmap='gray'); plt.axis('off'); plt.title('g, Shifted');
plt.subplot(1,2,2); plt.imshow(circle, cmap='gray'); plt.axis('off'); plt.title('g, Unshifted');
# Create an image with just a few non-zero pixels.
f = np.zeros([256,256])
f[50,100] = 1.0
f[48,110] = 2.
plt.figure(figsize=(8,8))
plt.imshow(f, cmap='gray'); plt.axis('off'); plt.title('f');
g = ifft2( fft2(f) * fft2(circle))
plt.figure(figsize=(8,8))
plt.imshow(abs(g), cmap='gray'); plt.axis('off');
plt.title(r'$(f \ast g)$');
f = plt.imread('pd.jpg')
f = f[:,:,0]
plt.imshow(f, cmap='gray'); plt.axis('off'); plt.title('An image f');
g = ifft2( fft2(f) * fft2(circle))
#plt.figure(figsize=(8,8))
plt.imshow(abs(g), cmap='gray'); plt.axis('off');
plt.title('g convolved with the image');
# Make and edge detector
g = Gaussian2D([256,256], 20)
G = ifftshift( fft2(fftshift(g)) )
ramp = np.outer(range(256), range(256)) - 128**2
H = G * ramp*1.j
circle = np.real( ifft2( fftshift(H) ) )
plt.figure(figsize=[8,5]);
plt.subplot(1,2,1); plt.imshow(fftshift(circle), cmap='gray'); plt.axis('off'); plt.title('Edge Kernel (Shifted)');
plt.subplot(1,2,2); plt.imshow(circle, cmap='gray'); plt.axis('off'); plt.title('Edge Kernel (Unshifted)');
f = plt.imread('pd.jpg')
f = f[:,:,0]
F = fftshift( fft2( f ) )
plt.figure(figsize=[10,20])
plt.subplot(1,2,1)
plt.imshow(np.log(abs(F)+1), cmap='gray'); plt.axis('off'); plt.title('Shifted DFT (log-modulus)');
plt.subplot(1,2,2)
plt.imshow(f, cmap='gray'); plt.axis('off');
g = fftshift( ifft2( fft2(ifftshift(f)) * fft2(circle)) )
plt.figure(figsize=(16,8))
plt.subplot(1,2,1);
G = fftshift(fft2(ifftshift(g)))
plt.imshow(np.abs(G), cmap='gray'); plt.axis('off')
plt.title('DFT(Edge Kernel) x DFT(Image)');
plt.subplot(1,2,2); plt.imshow(abs(g), cmap='gray'); plt.axis('off');
plt.title('Edge Kernel convolved with image');
thresh = 10
rows, cols = np.shape(f)
ctr = np.floor(np.array(np.shape(f))/2) +1
rr, cc = np.mgrid[-ctr[0]:(rows-ctr[0]),
-ctr[1]:(cols-ctr[1])]
def PlotFiltImg(F, filt):
F_filt = (filt)*F
f_filt = ifft2( fftshift(F_filt) )
plt.figure(figsize=[15,15])
plt.subplot(2,2,1); plt.imshow(f, cmap='gray'); plt.axis('off');
plt.subplot(2,2,3); plt.imshow(np.log(abs(F)+1), cmap='gray'); plt.axis('off');
plt.subplot(2,2,4); plt.imshow(np.log(abs(F_filt)+1), cmap='gray'); plt.axis('off');
plt.subplot(2,2,2); plt.imshow(np.real(f_filt), cmap='gray'); plt.axis('off');
filt_mask = abs(cc)<thresh
filt = np.zeros(np.shape(F))
filt[filt_mask] = 1.0
PlotFiltImg(F, filt)
filt_mask = abs(rr)<thresh
filt = np.zeros(np.shape(F))
filt[filt_mask] = 1.0
PlotFiltImg(F, filt)
filt_mask = np.sqrt(rr**2 + cc**2)<thresh
filt = np.zeros(np.shape(F))
filt[filt_mask] = 1.0
PlotFiltImg(F, filt)
Applying that filter in the frequency domain is the same as convolving the image with the kernel below, attained by taking the IDFT of the filter kernel.
Filt = fftshift(ifft2(ifftshift(filt)))
plt.figure()
plt.subplot(1,2,1); plt.imshow(filt, cmap='gray'); plt.title('DFT of filter'); plt.axis('off')
plt.subplot(1,2,2); plt.imshow(np.log(abs(Filt)+1), cmap='gray'); plt.axis('off'); plt.title('Spatial filter');
You can scale an image using the DFT by padding in cropping in the two domains (time/space, and frequency).
# We will magnify an image by padding its DFT coefficients, and then cropping the image
# that results from the IDFT of those padded coefs.
scale_factor = 1.5
padamt = int(rows*(scale_factor-1.)/2.)
G = np.pad(fftshift(fft2(f)), (padamt,), mode='constant')
g512 = ifft2( ifftshift(G) )
if padamt!=0:
g = g512[padamt:-padamt,padamt:-padamt]
else:
g = g512
plt.figure(1,figsize=[15,15])
plt.clf()
plt.subplot(1,2,1)
plt.imshow(f, cmap='gray'); plt.title('Original')
plt.subplot(1,2,2)
plt.imshow(np.real(g), cmap='gray'); plt.title('Scaled x'+str(scale_factor));
You can shift an image by adding a linear 'ramp' to the phase portion of its DFT coefficients.
dr = 20
dc = 12
ramp = -2.*np.pi*1j*(cc*dc + rr*dr)/rows
ramp = np.roll(np.roll(ramp,-1, axis=0),-1, axis=1)
wave = np.exp(ramp)
plt.figure()
plt.imshow(np.imag(ramp))
plt.title('Ramp in imaginary part');
# Make sure the centre pixel of the ramp has a phase of zero.
ctr = int(np.floor(256/2))
print('This should be zero -> '+str(ramp[ctr,ctr]))
print('This is the DC of our image -> '+str(F[ctr,ctr])) # This should be the DC of our image.
This should be zero -> -0j This is the DC of our image -> (2491737+0j)
G = wave * fftshift(fft2(f))
g = ifft2(ifftshift(G))
plt.figure(1,figsize=[15,15])
plt.subplot(1,2,1)
plt.imshow(f, cmap='gray'); plt.title('Original')
plt.subplot(1,2,2)
plt.imshow(np.real(g), cmap='gray');
plt.title('Shifted (down,right)=('+str(dr)+','+str(dc)+')');
You can detect a shift between two similar images by looking for a ramp in their phase difference. Or, you can take that phase difference and take the IDFT of it, and look for a spike. The offset of the spike from the origin gives the relative shift.
plt.figure(figsize=(8,8));
g = ndimage.shift(f,[dr,dc]) + np.random.normal(0.0, 5, size=np.shape(f))
G = fft2(g)
F = fft2(f)
H = G/F
plt.imshow(np.real(fftshift(H)), cmap='gray', vmin=-4, vmax=4);
%matplotlib notebook
plt.figure(figsize=(8,8));
h = ifft2(H)
plt.imshow(np.real(h), cmap='gray');
# Where does that lonely pixel occur?
coords = np.unravel_index(np.argmax(np.real(h)), h.shape)
print('Bright pixel at '+str(coords))
print('True shift was '+str((dr,dc)))
Bright pixel at (20, 12) True shift was (20, 12)
You can also rotate an image by simply applying the same rotation to the frequency domain.
theta = 30
f = plt.imread('pd.jpg')
f = np.array(f[:,:,0], dtype=float)
padamt = int(rows/2)
#padamt = 0
# Padding the image first makes the frequency domain smoother
ff = f
ff = np.pad(f, (padamt,), mode='constant')
FF = fftshift( fft2( ifftshift(ff) ) )
# Apply rotation to the Fourier coefficients
Gr = ndimage.rotate(np.real(FF), theta, reshape=False, order=1)
Gi = ndimage.rotate(np.imag(FF), theta, reshape=False, order=1)
G = Gr + 1j*Gi
g = fftshift( ifft2( ifftshift(G) ) )
if padamt!=0:
g = g[padamt:-padamt,padamt:-padamt]
%matplotlib inline
plt.figure(figsize=[15,15])
plt.subplot(1,2,1)
plt.imshow(f, cmap='gray'); plt.title('Original'); plt.axis('off')
plt.subplot(1,2,2)
plt.imshow(np.real(g), cmap='gray');
plt.title('Rotated by '+str(theta)+' degrees'); plt.axis('off');
The Fourier projection theorem is the basis for reconstruction in Computed Tomography (CT) scans.
f = plt.imread('ct.jpg')
f1 = np.array(f[:,:,0], dtype=float)
plt.figure(1, figsize=(10,10))
plt.clf()
plt.subplot(2,2,1)
plt.imshow(f1, cmap=plt.cm.gray);
px = np.sum(f1,axis=0)
plt.subplot(2,2,3)
plt.plot(px)
plt.axis('tight');
F = fftshift(fft2(f1))
plt.subplot(2,2,2)
Fnice = np.log(abs(F)+1)
plt.imshow(Fnice, cmap=plt.cm.gray);
ctr = GetCentre(F)
Fx = F[ctr[0],:]
plt.plot([1,254],[128,128], 'r--')
plt.axis([0,255,255,0])
plt.subplot(2,2,4)
plt.plot(np.log(abs(Fx)+1), 'r')
plt.axis('tight');
Px = fftshift( np.fft.fft(px) )
#plt.subplot(2,2,4)
plt.plot(np.log(abs(Px)+1), 'b--');
x, y = np.mgrid[0:256,0:256]
f = np.sin(2.*np.pi*x/256*40)
#f = ndimage.imread('pd.jpg')
plt.figure(figsize=(15,7))
plt.clf()
plt.subplot(1,3,1)
plt.imshow(f, cmap='gray');
c = plt.axis()
plt.title('Original');
g = f[::4,::4] # every fifth sample
plt.subplot(1,3,2)
plt.imshow(g, cmap='gray'); #plt.axis(c);
plt.title('Subsampled');
plt.subplot(1,3,3);
plt.imshow(g, cmap='gray'); plt.axis(c);
plt.title('Subsampled');
q = plt.imread('bricks.jpg')
q = np.array(q[:,:,0], dtype=float)
nr, nc = np.shape(q)
plt.figure(figsize=(15,7));
plt.subplot(1,3,1)
plt.imshow(q, cmap='gray', interpolation='none')
plt.title('Original'); c = plt.axis();
q_small = q[::5,::5]
qsr, qsc = np.shape(q_small)
plt.subplot(1,3,2)
plt.imshow(q_small, cmap='gray', interpolation='nearest')
plt.axis(c);
plt.title('Subsampled, with aliasing');
plt.subplot(1,3,3)
plt.imshow(q_small, cmap='gray', interpolation='nearest')
plt.title('Subsampled, with aliasing');
plt.figure(figsize=(15,7));
plt.subplot(1,3,1)
plt.imshow(q, cmap='gray', interpolation='none')
plt.title('Original');
q_small = q[::5,::5]
qsr, qsc = np.shape(q_small)
g = Gaussian2D(q.shape, 10) # Create filter
Q = fftshift( fft2(q) ) # DFT of the filter
G = fftshift( fft2(fftshift(g)) )
h = ifft2(ifftshift(Q*G)) # Apply filter in the freq domain
plt.subplot(1,3,2)
plt.imshow(np.real(h[::5,::5]), cmap='gray', interpolation='nearest')
plt.axis(c);
plt.title('Smoothed and subsampled (no aliasing)');
plt.subplot(1,3,3)
plt.imshow(np.real(h[::5,::5]), cmap='gray', interpolation='nearest')
plt.title('Smoothed and subsampled (no aliasing)');